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Abstract. This paper considers the extreme type-II Ginzburg-Landau equations that model 
vortex patterns in superconductors. The nonhnear PDEs are solved using Newton's method, and 
properties of the Jacobian operator are highlighted. Specifically, it is illustrated how the operator 
can be regularized using an appropriate phase condition. For a two-dimensional square sample, the 
numerical results are based on a finite-difference discretization with link variables that preserves the 
gauge invariance. For two exemplary sample sizes, a thorough bifurcation analysis is performed using 
the strength of the applied magnetic field as a bifurcation parameter and focusing on the symmetries 
of this system. The analysis gives new insight in the transitions between stable and unstable states, 
as well as the connections between stable solution branches. 
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1. Introduction. In this article, we study the symmetry-breaking transitions 
between stable and unstable patterns in small-sized superconducting samples. Su- 
perconductors are materials that expel magnetic fields and exhibit zero electrical 
resistance when they are below a characteristic temperature T^. Mathematically, 
the superconductor's states are described by a set of nonlinear PDEs, known as the 
Ginzburg-Landau system [21]. 

For simplicity, let us suppose that a sample of superconducting material occupies 
an open, bounded region of the EucHdean space, immersed in an external magnetic 
field Ho (see Figure 1.1). Above a critical temperature Tc, the material behaves like a 
normal conductor: it exhibits electrical resistivity and is homogeneously penetrated 
by the applied magnetic field. The material is said to be in a homogeneously non- 
superconducting state (or normal state). 

At low temperatures, T < Tc, the material exhibits a complete loss of resistivity, 
resulting in the formation of superconducting currents in the sample. Such currents 
give rise to an induced magnetic field and the total magnetic field B is expelled from 
the interior of the sample. Below a certain critical field strength i^ci, the magnetic 
field is expelled entirely; the material is said to be in a homogeneously superconducting 
state. Below the critical temperature and for stronger applied magnetic fields, however, 
mixed configurations can exist: the magnetic field penetrates only in confined regions 
of the sample. For so-called type-II superconductors [25], those areas are circular 
vortices, arranged in characteristic patterns. 

In large samples, the vortices organize in a regular pattern, also known as the 
Abrikosov lattice (see [1], [25] and references therein). In small samples, however, owing 
to the boundaries, the observed patterns can significantly deviate from the regular 
lattice and their organization depends sensitively on the intensity of the applied 
magnetic field as well as the geometry and the symmetries of the sample. These 
small-scale (mesoscopic) systems with simple geometric shapes like discs, triangles, or 
squares are of technological interest since they can be built into nanoscale devices [5]. 
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Fig. 1.1: States of a superconducting sample immersed in an external magnetic field Hq. 
Top: below the critical temperature and for H < i^ci, the sample is in a homogeneously 
superconducting state in which internal currents are generated and the total magnetic 
field B is expelled from the specimen (left); above the critical temperature or for 
H > i^c2, the material is in a normal state and the external magnetic field penetrates 
the whole sample (right). Bottom: type-II superconductors can exhibit mixed states, in 
which vortices of normal conductivity are embedded in a superconductive background. 
In the mixed configuration, B can penetrate the sample only through the vortices, 
giving rise to characteristic superconductive patterns. 



In applications, one is interested in finding steady states of the system, studying 
their stability and their dependence upon the external magnetic field. The state of 
a superconducting sample is, in general, characterized by two quantities: the total 
magnetic field B = : ^ and the density pc • ^ U dft ^ R of electron pairs which 
constitute superconductivity (Cooper pairs). 

A typical approach for studying superconducting states is to define a suitable 
Gibbs energy for the system and to derive a set of evolution equations for the order 
parameter i/j: QUdQ. C, IV^p = pc, and the magnetic vector potential A: R^ ^ R^, 
V X A = B. The resulting system in known as the Ginzburg-Landau system [21]. 
The associated initial-boundary- value problem has been studied both analytically an 
numerically. Various results on the existence and uniqueness of solutions, for example, 
can be found in [10, 34, 32] and references therein. 

However, it is often necessary to resort to numerical simulation to study the 
complex interaction of vortices in samples of arbitrary shapes: a popular strategy 
is to time-step the Ginzburg-Landau system via Gauss-Seidel iterations until an 
equilibrium is reached; the external magnetic field is then varied quasi-statically, and a 

2 



-1 

0.5 1 1.5 2 0.5 1 1.5 2 0.5 1 1.5 2 

(a) Disk. (b) Square. (c) Triangle. 

Fig. 1.2: (Reproduced from [9].) Typical cascades of branches of steady states of the 
Ginzburg-Landau problem for various two-dimensional sample shapes. The plots show 
the strength of the applied magnetic field (homogeneous and perpendicular to the 
sample) versus the normalized energy of the states. 



new steady state is found [9, 35]. We show a typical result of this analysis in Figure 1.2. 
The solution branches appear disconnected: when an instability is met, the direct 
simulation jumps to a nearby stable branch, as the employed numerical method can 
compute only stable solutions. 

The plots in Figure 1.2 are in good agreement with the hysteretic behavior that 
has also been observed experimentally [36], but they are not yet fully understood 
from the point of view of bifurcation analysis. The main results in this direction are 
confined to one-dimensional spatial domains (see [17, 4, 2] and references therein). 
The two-dimensional case has been studied by means of direct numerical simulation 
for various material parameters and strengths of the applied magnetic field [3], as well 
as for a variety of different shapes and domain sizes (see [19, 9, 13] and references 
therein), but the bifurcation scenario of the Ginzburg-Landau problem in two and 
three dimensions is largely unexplored. 

The main motivation of the present paper is to classify the instabilities occurring in 
superconducting samples using numerical continuation, as opposed to time-dependent 
simulations. To this end, we define a well-posed boundary- value problem, choose a 
spatial discretization, and find steady states of the Ginzburg-Landau problem by 
Newton iterations. More specifically, we focus on square samples of extreme type-II 
superconductors subject to a homogeneous external magnetic field. In this case, the 
Ginzburg-Landau problem simplifies considerably as it is possible to derive the vector 
potential A explicitly and then solve a nonlinear partial differential equation for the 
order parameter ip. 

We expect that the symmetries of the problem will influence the bifurcation 
landscape. As we will see, the relevant groups for the computations presented in this 
paper are the circle group and the dihedral group 1^4. The discrete symmetry 
suggests that we can use the equivariant branching lemma to predict symmetries of the 
emerging branches at bifurcation points [24, 27]. On the other hand, the continuous 
S*^ -symmetry induces the presence of a zero eigenvalue in the spectrum of the linear 
operator associated with the boundary-value formulation, causing problems to the 
convergence of the Newton iterations. 

We regularize the system by extending the boundary- value problem and employing 
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a suitable phase condition, using the framework proposed by Champneys and Sandstede 
[14]. The extended boundary- value formulation is then discretized using a common 
gauge-preserving technique and the patterns are path-followed in parameter space via 
pseudo-arclength continuation. 

To the best of our knowledge, this approach has never been employed before 
for the Ginzburg-Landau problem, albeit the application of Newton's method has 
been proposed in [21] and inexact Newton methods are often used in practice [28]. In 
addition, equivariant bifurcation theory has never been used to explain the instabilities 
found experimentally and numerically in superconducting samples, even though the 
importance of symmetries was pointed out in [15], where the system is linearized 
around the trivial steady state and the relative eigenmodes are studied in the context 
of C4-symmetries. 

The main result of the present paper is a classification of the bifurcations occurring 
in square domains of small and moderate sizes. In small samples, where the domain 
can host just a single superconducting vortex, the bifurcations are entirely determined 
by the natural two-dimensional irreducible representation of D4 (see [27], Section 4.3). 
However, as the domain size increases, the bifurcation diagram gets more complicated 
and it involves also one-dimensional irreducible representations of D4. Furthermore, 
in larger samples we compute stable vortices of higher multiplicity: these structures 
were previously found by direct numerical simulation [9], but their formation was still 
an open problem; our analysis shows that vortices with different multiplicity are all 
linked in parameter space via symmetry-breaking bifurcations. Furthermore, we used 
Newton-Krylov methods to solve the system, exploiting the properties of the Jacobian 
operator in the Krylov iterations. 

The remainder of the article is organized as follows. Section 2 discusses the 
Ginzburg-Landau system in the large- />: limit and details its symmetries. Section 3 
contains material on the linearization of the Ginzburg-Landau system and its self- 
adjoint ness, which is of importance for the numerical solution of the linear system 
associated with each Newton iteration. Section 4 is concerned with the regularization 
of the equations. Details on the discretization of the system with link variables and 
properties thereof can be found in Section 5. The numerical computations are included 
in Section 6, where we show families of solutions as a function of the strength of the 
applied magnetic field. We relate the bifurcations to the symmetries of the solutions 
with the help of the equivariant branching lemma. The appendix contains an extension 
of Keller's bordering lemma which is used in Sections 4 and 5. 

Notations. Throughout this article we will use bold symbols (x) for vector-valued 
quantities. For any z G C, and ^{z) denote its real and imaginary parts, z is 
used for complex conjugation. Similarly, for t/j : Q ^ ip : Q ^ C is such that 
ip^x) := ip{x) for all x G O. For function spaces, we use C^{Q) denote the vector 
space of all k times different iable functions. We use L'^{^) for the Hilbert-space in 
the field K of square-integrable functions over il, equipped with the inner product 
((^, i/j) = Ip^l) dO. for all '0 G The range of a Hnear operator L is denoted by 

1Z{L). For the symmetry groups under consideration, the symbol is used to denote 
the circle group {z ^ ^ : \z\ = 1} (which is group-isomorphic to S0(2)). For a given 
state ?/^, denotes the symmetry group under the action of which is invariant. 

2. The Ginzburg-Landau equation. For an open, bounded domain (7 C M^, 
with a piecewise smooth boundary dVL^ the Ginzburg-Landau problem is usually 
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derived by minimizing the Gibbs free energy functional 



(2.1) 



where the state ('?/^, A) is in the natural energy space such that the integral is well- 
defined [21]. As we have seen in the introduction, the scalar ip is commonly referred 
to as the order parameter, while A is the magnetic vector potential corresponding to 
the total magnetic field. The physical observables associated with the state (?/^, A) are 
the density pc = of the superconducting charge carriers and the total magnetic 
field B = V X A. The constant Gn represents the energy associated with the entirely 
normal (non-superconducting) state. 

The energy (2.1) is written in its dimensionless form and it depends upon the 
impinging magnetic field Hq and the material parameters a, /3, ^ G M. The most 
relevant parameters are n and ^; in particular, = A/^ is the ratio of the penetration 
depth A (the length scale at which the magnetic field penetrates the sample) to the 
coherence length ^ (the characteristic spatial scale of ?/^). A superconductor is said to 
be of type I if k < l/\/2, and of type II otherwise. 

To complete our description of the Gibbs energy, we remark that we have scaled 
the domain Q in units of the coherence length ^ while another common choice is to 
scale the domain by A [21]. 

Starting from the Gibbs energy and using standard calculus of variations, it is 
possible to derive the Ginzburg-Landau equations [21], a boundary- value problem 
in the unknowns and A. As anticipated in the introduction, we will simplify the 
problem and consider only the limit ^ oo {extreme type-II superconductors): this 
approximation gives satisfactory results for all high-temperature superconductors with 
large but finite values of a^, typically 50 < < 100. 

In this case, the Ginzburg-Landau problem decouples and we have 



Since for this decoupled system there are no magnetization effects, the magnetic fields 
Ho and B coincide. 

Since the sample's width scales with ^, the Isnge-hi limit A ^ means that the 
magnetic field A penetrates the whole sample, independently of V^. 

In passing, we note that Equation (2.2) does not coincide with the so-called 
Complex Ginzburg-Landau equation (see [6] and references therein). 

In the present paper, we consider a two-dimensional square sample 




in 



(2.2) 



where A = A(Ho) is given by the relations 




(2.3) 



n = nd-= {{x,y,z) e : {x,y) e {-d/2,d/2f,z 



0}, deR+, 
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subject to a perpendicular, homogeneous magnetic field Hq = (0,0, /i)^, /i G M. From 
(2.3) we can derive an expression for the induced vector potential 

A(x,?/;/i) = {Ax{x,y; /j.),Ay{x,y] /j.))^ := ^(-/i?/, /ix)^, (2.4) 

where we have deliberately omitted the third component. 

In conclusion, we will consider the following boundary- value problem with 
being the natural energy space over associated with the Gibbs energy (2.1) and 
its dual space. The equations are 

gC{ilj;/j.) := < 

[n • (-iV - A(/i))?/^ on dQd, 

where A(/i) is given by (2.4), with the parameters /i G M, (i G and with the boundary 
conditions given in the sense of traces. To shorten the notation, the dependence of A 
on jii will often not made explicit in the remainder of the text. 

Note that, because fid is convex and A G (7^(1^), solutions in the natural energy 
space immediately have higher regularity [8] and in fact coincide with the classical 
strong solutions in C'^{Q) H C^{Ti). 

Symmetries. As mentioned in Section 1, symmetries play an important role in 
the bifurcations scenario of our problem. The Ginzburg-Landau system for extreme 
type-II superconductors, (2.5), is left invariant by the action of the circle group S^^ 

l9^:V^i — ^V^exp(i7^), 7^G[0,27r). (2.6) 

The circle-group symmetry is also referred to as phase symmetry. 
In addition, 6C(V^; /i) is invariant under rotations by 7r/2 

p: V^(x, y) I — > i^{-y, x) (2.7) 
and conjugated mirroring along the ^-axis. 



a: ipix.y) i — > ipi-x.y). (2.8) 

Note that, up to conjugation in a, these are the classical group actions that generate 
the D4 symmetry group of the square. In fact, the group generated by p and a is 
isomorphic to D4. Even though symmetries of the Ginzburg-Landau problem have 
been considered before [15], the analysis was limited only to a linearization of the 
Ginzburg-Landau operator in the presence of rotations (2.7); in our case, we will 
consider the nonlinear problem and account also for conjugate reflections (2.8). 

In conclusion, the relevant symmetry group for our problem is generated by the 
actions (2.6)-(2.8), 

r:={Or,,p,a)^S^ xD^. 

We refer to the reader to Section 4, where we will explain how to factor out the 
continuous 5^ -symmetry that induces a singularity in the boundary- value problem 
associated with gC{ip; p)^ and we conclude this section by showing in Figure 2.1 a few 
examples of patterns computed via numerical continuation. The Ginzburg-Landau 
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Fig. 2.1: Several families of stable solution branches of the Ginzburg-Landau equation 
(2.5) as a function of the bifurcation parameter /i, the intensity of the applied magnetic 
field. The results, obtained via numerical continuation, are computed for a square of 
side length d = b.b. Top: On the vertical axis, we plot a measure of the Gibbs energy 
of the states (see remark 2 on page 17). We observe a cascade of instabilities similar 
to the one shown in the center panel of Figure 1.2, albeit the number of solution 
branches is larger in the latter diagram. As we will see in Section 6.2, the number of 
stable primary branches increases with the domain size d (the computations in the 
center panel of Figure 1.2 are for d = 7.1). Bottom: Selected stable patterns along 
the branches. Vortices are characterized by a localized region of low supercurrent 
density (blue in the pictures), as well as a 27r/c-phase change in arg?/^, where k is the 
multiplicity of the vortex. For higher values of /i, vortices with higher multiplicities 
are found: pattern 3, for instance, has a 2 x 27r phase change, it is a giant vortex with 
multiplicity 2. All patterns shown have full D4-symmetry and we expect to interpret 
the diagram in terms of symmetry-breaking bifurcations. 
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problem (2.5) possesses two trivial solutions: 6£(0; /i) = for all /i G M (the normal 
state) and 6/2(1; 0) = (the homogeneously superconducting state). As expected, we 
find branches of nontrivial stable D4-symmetric solutions arranged in a characteristic 
cascade, in agreement with the results obtained by direct numerical simulations where 
the parameter fi is varied quasi-statically (see Figure 1.2). 

3. The Jacobian operator. The patterns shown in Figure 2.1 were computed 
as regular zeros of a nonlinear system of equations derived from the Ginzburg-Landau 
problem (2.5). The solutions were found via Newton-Krylov iterations, that require 
the specification of the action of the Jacobian associated with 6C(V^; /i) (see [30] for 
details on iterative linear solvers). Even though there were previous attempts to solve 
(2.5) with a modified Newton's method [28], those implementations did not retain 
second-order convergence. Before deriving explicitly the regularization procedure that 
allowed us to compute the superconducting patterns, we introduce in this section the 
Jacobian operator J7('0; /i) associated with 6C(?/^;/i), and prove its self-adjointness 
with respect to a suitably-defined inner product in X^. 

For a given d G M+, /i G M, and i/j^SiIj G X^, let us consider 

= (-iV - Af{tlj + 6iIj) - (V^ + 6tlj) (l - (V^ + (5^)(V^ + (5^)) 

- [(-iV - A)^il; - (1 - i^t/j)] 
= (-iV - AfSip {ip Sip ^ ipSip ^Sip Sip) 

-Sip{i- i^ip) 

-\- Sip [ip Sip -\- ip Sip -\- Sip Sip) . 
Neglecting higher-order terms in Sip, we obtain the Jacobian operator 

J{t. m)^ := ((-iV - A)2 - 1 + 2|V^|2) + 

Note that J{ip; fi) is indeed linear when defined over and as R- vector spaces. 

We are now going to prove that the Jacobian operator (3.1) is self-adjoint with 
respect to the inner product 

This property allows us to employ standard methods for symmetric hnear systems such 
as the conjugate gradient or the minimal residual method (using this inner product) 
to invert the Jacobian at each Newton iteration. Note that (•, coincides with the 
natural inner product in (L^(0))^, which is isomorphic to L^(0), because for any 
given pair (p^ip ^ L^(^7), one has 



^(P\ mip 



The following lemma gives insight into the adjoint of a linear operator of the 
form (3.1). The lemma is formulated for general Hilbert spaces, and we will use it 
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for H = L^(r^); in our case, the operation C mentioned below will be the pointwise 
complex conjugation. 

Lemma 3.1. Let H be a Hilbert-space with inner product (•, and let there be 
an operation C : H ^ H such that 

C{ax) = aC{x) for all a eR,x e H (3.3) 
3? {C{x),y)^ = 3? (x, C{y))H for all x.yeH. (3.4) 

Let Ci^ C2 H ^ H be linear operators. For every x e H, let 

Cx := Cix + C2C{x). 

Then C is a linear operator on H as M.-vector space, and its adjoint £* with respect to 
the inner product Or ^ ('5 ^'"^ given by 

where C{, li\ are the adjoint operators in H of Ci, C2, respectively. 
Proof. Let x^y G and consider 

(x, Cy)^ = 3? (x, Cry + C2C{y))H = 3? (x, Cry)^ + 3? (x, C2C{y))H . 

Using the operator adjoints £5, we get 

{x, Cy)^ = 3? {Ctx, y)H + 3? {C^x, C{y)) 

= SR {Clx, y)j, + SR {C{C*2x), y) = {Clx + C{Clx),y)^ = {C*x, y)^ . 

□ 

Lemma 3.2. Let A(/i) G C^„{D,), n G {2,3}. The kinetic energy operator 

JCifi) : Xd ^ Yd, 

'(-iV-A)V, inQd (3-5) 
n • (— iV — A)(p, on dVid 



/C(/i)(^ := 



is self-adjoint with respect to the inner product (•, ■) l2(^^^~^ over the subspace C X^ 
with n • (-iV — A)(/p = 0. 

Proof. This result immediately derives from the fact that 

/ V^(-iV - A)Vd^^ = / (-iV - A)V^(-iV - A)LpdVt-i [ ^n-{-iV-A)ip 

Jq Jq JdQ 

(3.6) 

for ah i/j e L'^iytd), ^ ^ see [7]. □ 

Corollary 3.3. For any t/j G Xd and A(/i) G C^n{^), the Jacobian operator 
J{iI)]Ijl) defined in (3.1) is linear and self-adjoint over X^ with respect to the inner 
product (3.2), (•, ■)^. 

Proof By lemma 3.2, the operator of Ji(V^;/i) := ((-iV - A)^ - 1 + 2|V;p) 
defined over X^ with respect to the L^(Q(^)-inner product is self-adjoint. It can easily 
be checked that the adjoint operator of J2{'^), defined by J2{i^)^ = ?/^^(/? for (/? G Xd, is 
given by J2*('0)^ — tjj'^if. Also note that the complex conjugation fulfills the conditions 
(3.3). Application of lemma 3.1 then states that the adjoint of J7('0; fi) '. Xd ^ Yd is 
given by 



J*(V^; ii)lp = J^if + = JW + V^V = J{t^ m)^ 
for ah (f G Xd, and thus J'*(V^; /i) = J7(V^; /i). □ 



4. Nullspace and regularization with a phase condition. As stated in 
the previous sections, our aim is to compute solutions to the Ginzburg-Landau 
problem (2.5) and continue them in the parameter ja. This can be done in principle 
by discretizing the Ginzburg-Landau operator and applying standard numerical 
continuation techniques. 

However, as we have seen in Section 2, the boundary- value problem /i) = 

is invariant under the actions of the group T = x D4, and continuous symmetries 
(such as the phase symmetry determined by the circle group S^) make the problem ill- 
posed. After discretization, this leads to numerical difficulties that make it principally 
impossible to compute accurate approximations to the original problem [18] (see 
Figure 5.1a). 

This problem is usually met in computations of relative equilibria, which are 
time-dependent solutions whose temporal evolution is governed by a symmetry of the 
underlying differential equations. Typical examples are traveling waves (translational 
symmetries) and spiral waves (rotational symmetry). The continuous symmetry 
induces a zero eigenvalue in the Jacobian associated with the boundary- value problem, 
and it is therefore not straightforward to use Newton's method to compute the desired 
pattern: each Newton iteration inverts the Jacobian evaluated at a given solution i/j 
and requires a regular linear operator. 

A generic strategy to compute relative equilibria and to remove the singularity is 
to extend the boundary-value problem by introducing an additional scalar unknown 
and closing the system by means of a suitably-defined phase condition [11, 14, 33, 12]. 
The new boundary- value problem is well-posed, and therefore Newton's method can 
find the solution and path- follow it as a function of the parameters. This regularization 
technique can be applied to the stationary patterns of the Ginzburg-Landau problem 
to factor out the action of the circle group S^. To the authors' knowledge, the removal 
of the singularity for the Ginzburg-Landau system has not been considered in literature 
before. The regularization adopted here is an application of the framework proposed 
in [14]. 

In order to regularize the Ginzburg-Landau problem, we look at the action of 
G alg(5'^), 7^ G M, on a state i/j. Note that the exponential map of is given by 
exp : alg(6'^) ^ 6*^, ^ 6>^, where is identified with the action 6>^?/^ = e^'^i/j on a 
state tjj. This yields 



^ryV^ = ^exp(^^^)V^ 



t=o 



^=0 



and indeed the function it/j^ is in the nullspace of the Jacobian for a solution {ip^^ fis) 
of (2.5): 



J{^,; /is)(i^s) = [(-iV - A)2 - 1 + 2\^,\^] {i^,) - i^l^, 

= (1 - IV^sl') (iV^s) - iV^s + - i^l'^, = 0. (4.1) 

It is then possible to amend the Ginzburg-Landau problem and factor out the 
action of 5*^. To this end, for fixed /i, we compute (V^,/?) as a regular zero of the 
extended operator 

(V^, 7?) ^ {gc{^', fi) - - v^o)) 
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where $ : ^ R is a suitable phase condition and ipQ a given reference state. In the 
Ginzburg-Landau setting, the natural choice is the functional 



$ : ^ R, 3?(i^o, ^ - ^o) (4.2) 

with a given reference state i/jq G L'^{Q.d)^ subject to mild conditions (see corollary 4.1). 
Hence, instead of (2.5), we will consider the extended problem 

, = QCMm,)--=[^^^^^^^^^^^^^^^^^y (4.3) 

Remark 1. The phase condition featuring in (4.3) is also a necessary condition 

for 

minllV^o-V^e^^ll^,^^^^. 

This selects, out of all physically equivalent candidate solution states i/je^^, those two 
which are closest and furthest from tpo in the L'^{Q(i)-'^orm. 

If tps G Xd is a solution of the original equations (2.5), then (exp(ixs)'0s7 0)^ 
with Xs •= ~ arg((?/;o, '0s)Lg(Qd)) is a solution of (4.3) as well. The Jacobian operator 
corresponding to the extended problem (4.3) is 

Jp{^,r];fi):XdXR^YdXR, 

JpV^^mN 11 = 1 - I- 



^) \ ^((V^0,^)Lg(J7,)) 

We expect that the dimension of the nullspace of the extended Jacobian (4.4) is 
lower than the one of J{il)] ij). This is guaranteed by Keller's bordering lemma [29] 
if d\mk.eT J {^jj] ji) = 1. However, in the case of the Ginzburg-Landau operator we 
will encounter degeneracies of higher order. In the appendix, we present a bordering 
lemma that can be applied in such cases (Lemma A.l) and that is used in the proof of 
the following corollary. 

Corollary 4.1. Let (i/js^/j^s) ^ X^ x R be a solution of the original Ginzburg- 
Landau equations (2.5) with ifjs ^ 0^ and let i/jq G L'^{Qd) such that (i^o^i^s) 0- 
Then 

dimkeiJpitlJs'.liis) < dimkei J {i/js] j^is)- 



Proof The corollary is a direct consequence of lemma A.l (page 27) so it suffices 
here to verify that it can be applied on ^p(?/^./i). First note that the phase condition 
^((?/^o, V^)l2(j^^)) is a linear functional over the R-vector space X^,. Furthermore, we 

have by (4.1) that span{i?/^s} ^ ker J7('0s; Ms)- Evaluating the phase condition at i?/^s 
yields 

by assumption. Moreover, corollary 3.3 states that J{i^s] Ms) = yJ^ii^s] Ms) with respect 
to the inner product (3.2). From this, it follows that 

7^( J(V^s; Ms)) = ker( J*(V^s; Ms))^ = ker(J(V^s; Ms))^, 
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so to show that b = —itps ^ ^{J) is suffices to show that b is not orthogonal to all 
of ker(J'(?/^s; Ms)) with respect to the inner product (3.2). This holds true for ?/^s 7^ 
since 

Thus, all conditions of lemma A.l are fulfilled and its application to /i) concludes 

the proof. □ 

In the remainder of this section, we show that the extended operator retains the 
symmetries of the Ginzburg-Landau system. 

Symmetries of the extended system. Given 7 G F, the symmetries 7 of the extended 
system are defined to act 



Lemma 4.2. The extended system (4-3) is equivariant exactly under all actions in 
r that leave i/jq invariant^ i.e., T = H T. 

Proof. We have to proof equi variance only for the generators p, a G T. For a given 
7 G n F, it has to be shown that 

-f[gC{tlj) -ir]tlj]\ ( I \'- rr (-( I \\ ^^(tV^) " ^^(tV^) A 

-((^o,^),g(o.))j =^^^^(^'^) = ^^^(^(^'^))= V^((^o,7^>.g(o.))j' 

which holds obviously true for the first component, owing to the F-invariance of QC. 
As for the second component, we have 

After a suitable change of variables, this is equivalent to show that 

S^((^o,^>ig(o,)) = 3(^(7"¥o)^det^?^d!^^ G L|(i^) (4.6) 

holds exactly for all 7 G H F. 

Firstly, let us show this equivalence for the cyclic subgroup C4 -< S^^. Given 
7 G C4 (and thus 7?/^o = V^o, deti?^ = 1), equation (4.6) obviously holds true. On the 
other hand, let us assume that equation (4.6) is valid and let us take a sequence of 
Dirac-(5 functions centered at (xo,^o), i^^^^ = ^^[3.^ y^y Then 



= (V^o(^o, 1/0)) = /hn ^ {J^ ^^^^'^ = 



lim ^ 



^ 7" VoV'^'^ dfi^ = 9 ((7- Vo) {xo, yo)) 



Since this can be done for any (xo,^o) G 1^ U dft^ we have ^(V^o) G C4. The same 
result is obtained for 5R(V^o) by taking ip^^"^ = i^^l^,^ ^q). We conclude that 3 C4 is 
also necessary for (4.6) to hold. 

The very same arguments can be applied to the conjugate reflection a, noting that 
deti^cr = — 1, and that the action of a changes the sign of the expression. □ 
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The choice of tpo must hence be such that it eliminates the phase invariance 
(according to corollary 4.1), and that it preserves the other symmetries of the system 
(according to lemma 4.2). The first condition is equivalent to demanding ^ (ipo^ips) 
which indeed is a rather mild condition that will be fulfilled, for instance, by ipo = 1 
for most scenarios considered later. Note that it is also possible to update i/jq in each 
Newton step to the current guess V^*^^^: For a solution ?/;s, let i/j^ = V^*^^^ + e*^^^; we have 



which is guaranteed to be nonzero for sufficiently small e^^^ if 7^ 0. Note that 
intermediate Newton steps might not exactly preserve the symmetries of the system, but 
the symmetry breaking is weak in the sense that symmetry is preserved at convergence, 
and will not do harm [27]. 

5. The discretized system. An important property of the full Ginzburg- 
Landau equations is its gauge invariance^ a generalization of the phase symmetry 
(2.6) to the case where A is not fixed (see section 3.1 in [21]). While the reduced 
invariance with fixed A is preserved under all consistent pointwise discretizations of 
the Ginzburg-Landau equations [20, 22, 28], ordinary finite-difference discretizations 
lead to systems that are gauge invariant only up to 0(/i), where h is the grid spacing. 
It is thus customary to reformulate the equations using techniques from lattice gauge 
theory. For the convenience of the reader, the new system will be presented in this 
section. We also show that, for appropriate phase conditions, all the symmetries are 
preserved in the extended discretized system. 

5.1. Formulation with link variables. Let us consider the functions 
Ux{x, y) := exp ^-i j A3,(^, y) d^^ , 
Uy{x, y) := exp ^-i j Ay{x, v) dv^ , 
with arbitrary, fixed xq, ^0 ^ It can be checked easily that 

^ -U,{x, y)-^iU,i^) = (-iV - A)V - i(V • A)^. 

ue{x,y} 

The Ginzburg-Landau equations (2.2) can then be written as 

= - U^£2{U,i^)-^i^{l-H) inf^, 

= n • on on. 

Uu only appears in the product Uuip which guarantees preservation of full gauge 
invariance [21]. 

5.1.1. Discretization. Let N > 0, for simplicity even, h = d/N, and let Qh = 
{h-{ij) I - N/2 < i < N/2, -N/2 <j< N/2} be a uniform grid on Q^. Furthermore, 
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let us consider the discretization ?/;^^^ G = C^^^^^^^^^^^. The ordinary five-point 
discretization of (5.1) with centered finite differences on the boundaries is given by 







Vi,j G{-7V/2,...,7V/2} (5.2) 



where the finite-difference operator D^x is defined by 



( 



for i 
for - 
for i 



N/2, 
< i 
-N/2, 



2 ^ ^ ^ 2 ' 



h-^{-2(ui%+^,^^X^+2^'^]) 

and hkewise for Dyy'il)^^\ with unknown ij)^^^ G X/^, where 

(/7f)),±i,, =exp(-i4-;±n^.(^^,0)) +0(/.^^ 

with J^;^' {Ax{-, Vj)) ^ /J.'^' Aa,(^, ifj) d^, and hkewise for {U^^^)x,j±i. The order p of 
the approximation depends on the quadrature method in use. It is easy to verify that 
the discretization (5.2) has order of consistency min{2,p — 2}. 
The discrete Jacobian operator J'^^'^ is defined similarly as 



(h) 



ih) 



• (h) 



ih) 



(5.3) 



and it is self-adjoint with respect to the scalar product This is a consequence 

of lemma 3.1 upon realizing that the operators D^x and Dyy are both self-adjoint with 
respect to the inner product (•, •)c in A consequence of this is that all 

eigenvalues of the operator (5.3) are real-valued. 

5.1.2. Symmetries of the discretized system. Now that the structure of the 
discretized operator is described, we review the symmetries of the associated boundary- 
value problem. Many of the results in this section can be borrowed from Section 4 on 
the continuous problem with little modification. For example, the discretized system 
(5.2) is invariant under 



pointwise for each 77 G [0, 27r). Let further the discrete symmetry operators p and a be 
defined by 



ih) 



(5.4) 
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It can easily be shown that the discretized system (5.2) is invariant under these actions. 

Just Hke in the continuous case (4.1), the discrete phase invariance induces a 
nontrivial nullspace of J^^^: 



j 



ih) 

j 



0. 



•2i 



(h) 



4^ 



(5.5) 



and hence j G kei J'^^\ips ;/is)- This makes it impossible to treat the system (5.2) 
with generic linear solvers at a solution V^s^\ as the system is then exactly singular. 
In the neighborhood of V^s^\ the system will have a large condition number, making 
round-off errors dominate the update term [18] which in turn flaws the next Newton 
step. This phenomenon is illustrated in Figure 5.1a. 

As suggested in section 4, we will avoid the singularity of the Jacobian by using a 
phase condition, which in its discretized form reads 




N/2 

(^0 

= -N/2 



ih) 



,4';' 



(5.6) 



where ipQ^-' G is a given reference state, and 



Wi 



for -N/2<i< N/2, 
for i e {-N/2, N/2}, 



The discretized version of the extended system is then 




(5.7) 



and the symmetry operations for this extended system can be defined just like in (4.5). 
Parallel to (4.4), the discrete extended Jacobian operator j'-'^^ is 



J^('^)(V('^);r/;M) 



(h) 



/( ^) _ iry)<^W - (5.8 
J('^)(vr,<p('^)) 



The following two statements are the discrete versions of the lemmas 4.1 and 4.2. 
Corollary 5.1. Let g x R 6e a solution of the original dis- 

cretized Ginzhurg-Landau equations (5.2) with ipi^^ ^ 0^ and let ipQ^^ G X^ such that 
^0. Then 
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Fig. 5.1: The residual of the Newton iterations and the condition number hci of the 
associated linear system in the 1-norm. (a) Here, the linear Jacobian systems are solved 
using Gaussian elimination and deliver flawed Newton updates as the condition number 
increases. This is a principle problem and it is not limited to Gaussian elimination, 
(b) The regularization removes the singularity and this leads to bounded condition 
numbers. The linear systems can then be solved accurately. 



Proof. The proof runs parallel to the one of corollary 4.1, using the discrete inner 
product 5R(-, -)c' □ 

Lemma 5.2. The extended equations (5.7) are equivariant exactly under S -(h) H 

Proof. Again, the proof is essentially parallel to the one of lemma 4.2; instead of 
series of Dirac-J function (^(3^0,^/0)' their discrete equivalents 

^^'-^ 1^0 otherwise J ^'-^ 

□ 

6. Numerical results. Using the framework presented in the previous sections, 
it is possible to solve the Ginzburg-Landau equations numerically for any given 
parameter /i (the strength of the apphed magnetic field) and d (the edge length of the 
sample). As discussed in Section 1, the intensity of the applied magnetic field, can be 
tuned experimentally, and it is thus interesting to explore the bifurcation scenario as 
this parameter is varied. 

Because of the symmetries of the Ginzburg-Landau system posed on the square, 
we expect symmetry-breaking bifurcations to arise. As described in Section 4, the 
extended system (4.3) does not bear the continuous iS^-symmetry such that the relevant 
symmetry group for our computations is D4. 

Symmetry-breaking bifurcations in D4 are well known [27, 24]. We recall here 
that, in our case, the group generators are p, the rotation by 7r/2 (see equation (2.7)), 
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and cr, the conjugated mirroring along the ?/-axis (see equation (2.8)). We expect 
that symmetry-breaking bifurcations will occur when critical eigenvalues become 
unstable with (algebraic and geometric) multiplicity either 1 or 2. With a simple 
unstable eigenvalue, one should expect either a symmetry-preserving turning point 
or a pitchfork bifurcation with branches corresponding to the four one-dimensional 
irreducible representations of D4. With an eigenvalue of multiplicity 2 crossing the 
origin, two families of branches emerge from the bifurcation point, corresponding to 
the conjugacy classes of the D4 isotropy subgroups (p) and (crp), respectively. 

In the present section, the parameter ja will be varied for two different domain 
sizes d. The simplest nontrivial example of symmetry-breaking bifurcation occurs for 
small domain sizes, so we have deliberately chosen d = 3.0, a domain size that is just 
enough to host a single vortex. Subsequently, we study the case d = 5.5, for which the 
bifurcation scenario becomes increasingly more involved. 

The bifurcation diagrams are traced via standard numerical continuation methods 
[31]. The technical implementation is based on the Trilinos project [26] and exploits 
the sparse structure of the discrete Jacobian operator (5.3) as well as its properties, as 
outlined in Section 3. 

In the remainder of this section, we will denote solution branches (and relative 
patterns) alphabetically and bifurcation points with numerals. 

Remark 2. Unless otherwise stated, the bifurcation diagrams are plotted in terms 
of the expression 



(6.1) 



which is part of the Gibbs energy (2.1). This is in accordance to what is usually done 
in the physics literature. Applying (3.6) in the case t/j G X^, we obtain 



[ + [ ^/'(-iVV-A(^))2^dn . 

Only solutions ip{^) of the Ginzhurg-Landau equations (2.2) are considered, so that 



F{i,,n)=F-lJ-^ 



Thus, computing the significant portion (6.1) of the Gibbs energy (2.1) effectively 
reduces to evaluating 



6.1. Small-sized system {d = 3). The first computed solution corresponds to 
a superconductor in the absence of a magnetic field, that is, /i = 0; the system is in 
the homogeneous solution = 1 and it is said to be in a completely superconducting 
state. The solution has all the symmetries of the system (its isotropy subgroup is the 
full group D4) and is stable as a global minimum of the free energy (2.1). 

With the help of numerical continuation, a series of solutions for increasing /i is 
constructed. This results in branch A in Figures 6.1 and 6.2, showing the energy of 
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Fig. 6.1: Free energy of the solutions as a function of the strength of the apphed 
magnetic field /i. Solid (dashed) Hnes represent stable (unstable) states. Branch A with 
D4 symmetry starts at the homogeneous state with zero field and becomes unstable 
at bifurcation point 1. Branches G and C emerge from the bifurcation point, both 
characterized by a single vortex entering the domain; branch C has mirror symmetry 
along either the horizontal or the vertical center line, while branch G has mirror 
symmetry along one of the diagonals of the square domain. At bifurcation point 6, 
branches G and C connect to F, characterized by solutions with a single vortex in the 
center of the domain and full i^4-symmetry. 
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Fig. 6.2: The two largest eigenvalues of the Jacobian as a function of the applied 
magnetic field. We plot the eigenvalues close to bifurcation points 1 and 6, for each of 
the four solution branches in Figure 6.1. For field strength above fi = 1.646 the main 
branch is unstable, while for fields strengths weaker than /i = 1.175 the branch with a 
single vortex is unstable. We see that the largest eigenvalue of the main branch, which 
has multiplicity two, splits into two separate eigenvalues. Curve C has two unstable 
eigenmodes, while B has one stable and one unstable eigenmode. The colors reflect 
the branches in Figure 6.1. 



the solution and the two most unstable eigenvalues of the Jacobian as a function of 
the field strength, respectively. For non-zero field strength, the solutions deviate from 
the homogeneous superconducting state, developing zones of low supercurrent density 
near the edges of the domain (see pattern A in Figure 6.1). As ja is increased, the 
states are characterized by a higher energy and they maintain full D4 symmetry. 

At field strength fi ^ 1.64 (point 1 in Figure 6.1), an eigenvalue with multiplicity 2 
becomes unstable. At this bifurcation point, one can apply the equivariant branching 
lemma: the Ginzburg-Landau equation is equivariant under the symmetries of the finite 
group D4 and the eigenvalues cross the origin with non-zero speed, (see Figure 6.2). 
The lemma guarantees the existence of two solution branches emerging from the 
bifurcation, corresponding to the conjugacy classes of the isotropy subgroups (a) and 
(crp). They both have a one-dimensional fixed-point subspace. Hence, we expect two 
different families of solution branches, each containing four equivalent bifurcation 
curves with states belonging to one group orbit. The two families are found in the 
branches G and C of Figure 6.1. 

Before describing curves G and C, the two curves that emerge from the bifurcation 
point, we continue to follow the original branch A for increasing /i. The state is now 
unstable and retains full D4 symmetry. The magnetic field penetration increases from 
the boundaries until, at /i 1.89, the branch connects to the trivial state '0 = 0, which 
corresponds to the normal state of the sample. 

We now discuss the curves G and C, which have reduced symmetry and emerge 
from the bifurcation points 1 and 6. Curve C corresponds to solutions in which a 
single vortex moves in from one of the four sides of the square. These solutions belong 
to the conjugacy class of the subgroup (a) and the single vortex sits either on the 
horizontally or vertically centered line. 

The other family of solutions, on branch G, also features a single vortex entering 
the system, but along one of the diagonals. These solutions have an isotropy subgroup 
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Fig. 6.3: The four main branches found for d = 5.5. The corresponding stable patterns 
with one vortex in the middle of the domain are presented in Figure 2.1. Solid (dashed) 
lines represent stable (unstable) states. Shaded areas are detailed in Figures 6.5 and 
6.6. 



that belongs to the conjugacy class of (crp), hence their symmetry with respect to one 
of the diagonals. 

Solutions belonging to curves G and C are energetically similar, the latter having 
slightly higher energy, as it can be seen from the inset of Figure 6.1. 

As we decrease the field strength from point 1 to point 6, the vortex moves along 
the center line for curve C, or along the diagonal for curve G, towards the center of 
the sample. At field strength /i ?^ 1.18 (point 6 in Figure 6.1), each solution features 
a vortex in the middle of the domain and enjoys full D4-symmetry. As we can see 
in Figures 6.1 and 6.2, bifurcation point 6 is analogous to bifurcation point 1, but it 
involves branch F instead of A. 

The solution curve F in Figure 6.1 is characterized by a single vortex in the middle 
of the domain and is unstable for field strength weaker than /i ^ 1.18. This solution 
branch extends all the way up to field strength /i ^ 2.30 where it connects to the 
trivial zero solution. 

In a physical experiment where the magnetic field is first increased and then 
decreased, we would expect to observe hysteresis: while increasing, the system would 
initially follow branch A, switching to F at point 1; conversely, for decreasing /i, we 
would pass from branch F to A, at point 6. Hysteresis effects such as this one have 
been discussed in [3], and observed experimentally in many superconducting systems 
(see also Figure 1.2). 

6.2. Larger domain size (d = 5.5). In this section, we repeat the numerical 
experiment of Section 6.1 for a larger sample. In this context, it will be interesting to 
observe how the states of branch A destabilize: with edge length d = 5.5, more vortices 
can enter the domain, leading to a much more complicated bifurcation diagram. 
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Fig. 6.4: Schematic of solution branches, bifurcations, and patterns found for d = 5.5. 



Before starting to describe ah the branches found by means of numerical continua- 
tion, we anticipate that we found four main branches, as opposed to the case d = 3.0, 
where we computed only two. The four main branches are collected in Figure 6.3: 
they are labeled F, M, A, corresponding to states with vorticities 1, 2, 3 and 4, 
respectively. Their stable segments, together with a few corresponding patterns, have 
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previously been sketched in Figure 2.1. 

In the remainder of this section, we will concentrate on the two shaded areas (zone 

1 and II) of Figure 6.3. In these regions, a series of symmetry-breaking bifurcations 
and cross-connecting branches are found. 

As in the previous section, we start from the trivial homogeneous state tj; = 1 
at /i = 0, and increase /i. The resulting solution branch, enjoying full D4 symmetry, 
is labeled A and features four vortices entering the domain from the sides, similarly 
to what happens for d = 3.0. While this scenario resembles the one described in 
Section 6.1, the bifurcations occurring in zones I and II are quite different from the 
small-sized case, and we discuss them one by one in the remainder of this section. We 
refer the reader to the schematic in Figure 6.4, where we present all the branches, 
bifurcations, and representative patterns computed for d = 5.5. 

6.2.1. Zone I. Branch A in zone I destabilizes with a simple eigenvalue, at field 
strength ja ^ 0.70 (see point 1 in Figure 6.5). This mechanism is different from what we 
found the small-sized system, where an eigenvalue with multiplicity 2 becomes unstable. 
We can still apply the equivariant branching lemma: we expect a single family of 
solutions bifurcating from point 1, corresponding to a one-dimensional irreducible 
representation of D4 [27]. 

The corresponding branch is labeled B in Figure 6.5. It belongs to the conjugacy 
class of the isotropy subgroup {p^^cj)^ representing the mirror symmetries along 
horizontal and vertical center lines. When we follow this branch for decreasing values 
of /i, two vortices move simultaneously into the domain from opposite edges (left-right 
or top-bottom). 

Along branch 5, we find another symmetry-breaking bifurcation, point 2, where 
a second simple eigenvalue becomes unstable. This is shown in detail in the bottom 
panel of Figure 6.5, where we plot the negative value of the L2(r2/i)-norm in order to 
visualize the branches better. Branch C, emerging from point 2, has further reduced 
symmetry, corresponding to the conjugacy class (a), that is, a family of branches with 
a single vortex on one of the center lines, away from the center. 

On branch C, the vortex moves towards the middle of the sample and is connected 
via point 6, at /i ^ 0.25, to branch F, the second main branch with full 1^4 symmetry. 
A single vortex sits in the center of the domain throughout branch F and solutions on 
F are unstable for fields weaker than fi ^ 0.25. This branch is similar to branch F in 
the small system described in the previous section. 

Bifurcation point 6 features a null eigenvalue with multiplicity 2 and has the 
same symmetry properties as the bifurcation points discussed in Section 6.1. There, 
eigenvalues with multiplicity 2 became unstable on a branch with symmetry and 
two branches emerged with with symmetries {a) and {ap) (see also Figure 6.1). In the 
current system, it has already been found that branch C with symmetry (a) connects 
to point 6, and a second branch with symmetry (crp) is to be expected. This branch 
has a single vortex on one of the diagonals and is shown as curve G in Figure 6.5. In 
contrast to the small size system, this curve does not connect to bifurcation point 1. 
Instead, it connects to curve D via bifurcation point 5. 

A branch for which there is no equivalent in the smaller system is branch D in 
Figure 6.5, with a single vortex with multiplicity two (and hence phase change of 

2 X 27r, a so-called giant vortex)^ in the middle of the domain. Branch D has full 
D4 symmetry and is only stable for fields larger than p ^ 0.64. The corresponding 
bifurcation is marked by point 3 in Figure 6.5 and connects to branch B (see above). 

At point 3, the two vortices of B merge into the giant vortex; similarly, branch G, 
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Fig. 6.5: Bifurcation diagrams and representative patterns found in zone I of Figure 6.3, 
in the case d = 5.5. Top and middle panels: free energy versus magnetic field intensity. 
Bottom panel: the negative norm of the solution is used in the bifurcation diagram, 
in order to separate points 1 and 2. Solid (dashed) lines represent stable (unstable) 
states. 
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which emerges from point 6 on branch F, connects to branch point 5 on branch D. 

In the remaining part of zone I, we found that the main branch A has another 
instability, at bifurcation point 4. This bifurcation features a critical eigenvalue with 
multiplicity 2 and thus two families of solution branches emerge. Along branch 
three vortices enter the domain from three of the four sides of the domain. This 
branch corresponds to the conjugacy class of the subgroup (a). The three vortices 
move towards the center of the system along the branch where they finally merge into 
a giant vortex with multiplicity 3 at point 7, connecting to branch H. 

To conclude our exploration of zone I, we examined branch /, emerging from 
point 4 on the main branch A, for decreasing values of /i. Patterns on this branch have 
two vortices entering from two adjacent edges of the system. This branch is symmetric 
under reflections over one of the diagonals and corresponds to the conjugacy class of 
the subgroup (crp). 

6.2.2. Zone II. We now move to the upper part of the bifurcation diagram in 
Figure 6.1. An important difference from the small-sized system d = 3 is that the 
main branch A restabilizes as the fleld increases, as shown Figure 6.6. 

As we increase /i along the main branch A, four vortices are moving in from the 
midpoints of the edges towards the center; the solutions maintain full D4 symmetry. 
At field strength fi ^ 1.07, the four vortices arrive at the center and form a giant 
vortex with multiplicity 4. As the field strengthens further, this giant vortex breaks 
up again and four separate vortices move away from the center along the diagonals. 
Note that there is no bifurcation point associated with this reorganization as none of 
the eigenvalues of the Jacobian crosses the origin. 

At field strength /i 1.14, one of the unstable eigenvalues of bifurcation point 1 
restabilizes. This yields bifurcation point 8 in Figure 6.6. From 8, branch K emerges 
and connects to branch with a vortex of multiplicity 2 in the center of the domain. 
Along branch two of the four vortices are pushed out of the sample along one of 
the center lines, while the two remaining reorganize into a giant vortex of multiplicity 
2 (point 12). A sequence of patterns of branch K can be found in Figure 6.8. 

Branch A restabilizes at field strength /i ^ 1.15. The pattern with four symmetric 
vortices on the diagonals is now stable. The bifurcation point that marks this transition 
is labeled as point 9 in Figure 6.6. Two solution curves emerge from point 9, namely 
branches L and J. 

Branch L, connecting to branch M via point 10, features five vortices, as can be 
seen in Figure 6.9: four vortices arranged symmetrically, rather close to the center, 
and a single antivortex at the center of the domain, so that the total vorticity of the 
configuration is 3. A giant vortex of multiplicity 3 is formed at bifurcation point 10 on 
branch M, where it is unstable. The fact that the vortices do not arrange as a giant 
vortex with vorticity 3 in a stable fashion has been predicted in [15]. Solutions on M 
are unstable for weak fields strengths (see bifurcation 15 in Figure 6.6). 

In a similar way, branch J starts at point 9 and connects to point 11 on branch 
F for decreasing /i. The patterns along this branch are shown in the sequence of 
snapshots in Figure 6.7. 

At field strength /i ^ 1.50, the main branch A loses its stability again at point 13 
in a scenario similar to the small-sized system discussed in Section 6.1. The eigenvalues 
of the Jacobian at this bifurcation point are degenerate and two branches emerge, each 
of which has a single vortex entering either along the diagonals or along the center 
lines. These branches connect to a stable branch with five vortices organized like the 
five dots on a dice. This branch has been omitted in the figures. Further on the main 
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branch, a second simple eigenvalue becomes unstable at point 14. 

7. Discussion and conclusions. We have presented an initial exploration of 
the symmetry-breaking bifurcations of the vortex patterns as modeled by the Ginzburg- 
Landau equations. In the case of extreme type-II superconductors, we assumed a 
homogeneous applied magnetic field and showed how the vortices reorganize as the 
strength of the applied field is varied. In the small square domain = 3), we beheve 
to have given a complete account of the instabilities of the system. For a larger system, 
the bifurcation diagram becomes much more complicated, and we found a large number 
of states and symmetry-breaking bifurcations. 
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Fig. 6.7: Patterns of branch J in zone II (see Figure 6.6). Branch J bifurcates from 
branch F, which has a single vortex with multiplicity 1, and connects to branch A at 
bifurcation point 9. 



The paper also presents a study of the symmetries of the system. It has been 
shown that the continuous system bears symmetries isomorphic to S^xD4^. The 
discretization has been chosen in such a way that it preserves to machine accuracy 
both phase and geometric symmetries. 

Owing to the symmetries of the system, it is possible to use the Equivariant 
Branching Lemma in order to predict the existence of new branches at symmetry- 
breaking bifurcations, and subsequently compute them numerically. To the best of the 
authors' knowledge, most of the patterns contained in this paper are unknown to the 
physics community: even though unstable patterns can not be obtained experimentally, 
we point out that the methodology proposed in this context could be effectively used 
to find new stable patterns. 

The present paper analyzes the Ginzburg-Landau system on a square, but the 
same technique can be applied to all geometries with inherent symmetries, e.g., regular 
n-gons. It is not immediately obvious, though, how to choose the magnetic vector 
potential gauge such that the corresponding Ginzburg-Landau formulation remains 
equivariant with respect to Dn] some work in this area has been done in [16]. Note 
that, for increasing n, the ever more complicated subgroup structure of Dn will lead 
to different bifurcation scenarios [23, 24]. 

In the present paper we simplified the Ginzburg-Landau equations considering the 
large- Ai: limit, where the equation for the magnetic vector potential A decouples from 
the order parameter ?/^. It will be necessary, in the future, to study the bifurcations 
in the coupled system for intermediate and small values of n. However, this task 
will also pose new numerical challenges: the magnetic vector potential appears as an 
additional (vector- valued) unknown and its domain of definition is the whole space. 
In practice, the vector potential will approach its boundary condition defined by Hq 
sufficiently far away from the sample, but the validity of this approximation is still 
an open problem. The coupled system will in any case hold many more unknowns, 
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Fig. 6.8: Patterns on branch K in zone II (see Figure 6.6). The branch bifurcates off 
from branch A, with four vortices, at bifurcation point 8. Two of the four vortices are 
pushed out of the sample, while the remaining two reorganize into a giant vortex of 
multiplicity 2 at point 12. 



and a robust preconditioning strategy for solving the appearing Jacobian systems will 
be crucial. However, the regularization technique that we employed for the extreme 
type-II case is applicable for finite values of n and for generic spatial discretizations of 
the Ginzburg-Landau problem. 

Nevertheless, we believe that results of this paper are a first step in understanding 
the bifurcations in the coupled Ginzburg-Landau system for various mesoscopic systems 
that are relevant for nanoscale devices. The approach proposed here opens up the 
possibility of a systematic exploration of the solution landscape in regions that are 
precluded to direct numerical simulation. 
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Appendix A. Extension of Keller's bordering lemma. 

Keller's bordering lemma [29] provides conditions on how a finite-dimensional 
linear system with a singularity of dimension 1 can be regularized by adding an 
additional unknown as well as an additional equation. In the present context, however, 
it is necessary to formulate the lemma in general vector spaces. Also, the defect of 
the present problem may be greater than one. Such situations occur, for example, in 
several branch points described in section 6. The following lemma shows that it is 
always possible to remove one of the singularities. 

Lemma A.l. Let X, Y he 'K-vector spaces and let C : X ^ Y linear with 
dimker £ = k > 0. Let further b G Y , d G K, and f : X ^ K a linear functional. Let 
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Fig. 6.9: Patterns on branch L in zone II (see Figure 6.6). The patterns show five 
vortices: four vortices are arranged symmetrically, rather close to the center, and a 
single antivortex sits at the very center of the domain, so that the total vorticity of 
the configuration is 3. At bifurcation point 10, when a giant vortex of multiplicity 3 is 
formed, the pattern becomes unstable, on branch M . 



the operator £:XxIK^Y'xIK he defined by 



for all X = (x,^)^ G X X K. Then k := dimker>C < k if and only if b ^ ^[^) ^'^d 
there exists a v G kei C with f{v) ^ 0. 

Proof On the one hand, let b ^ 7^(>C) and let v G ker>C with f{v) 7^ 0. Let 
{{w'''^\^i)^}^^-^ C X xK denote a basis of ker>C, and take sl x e ker£. 



X — ^ a. 



i=l 

with arbitrary G K. With this representation, we have 

k k 



= /:^a,^(^)+6^a,e,, 

i=l i=l 

(k \ k 

i=l J i=l 

Because b ^ Tl{C), it must be Yli^i ^i^i = as otherwise 



-1 ~ 

k 



28 



Because the ai are arbitrary, we have = for ah i. Since {{w^'^^S^i)^}^^^ is 
hnearly independent in X x K and aU are zero, {w^'^'^}^^-^ is hnearly independent 

in X. Besides that, it follows that Yli=i^i^''^^ ^ keri2, and again because the ai 
are arbitrary, we have w^'^'^ G ker£ for all i G {1, . . . , Hence dimker£ > k. One 
can exclude k = dimker>C since then ker£ = span{^(;*^*^}^^^, and at the same time 
= f{Yli=i^i^^^'') arbitrary a^. This contradicts the assumption there is a 
V G ker£ with f{v) ^ 0. Hence k < k. 

On the other hand, let k < k. Consider the set W := keril x {0}. Obviously it is 
diml^ = A:, and additionally for any w = {w,0)^ G one has 

fC{w)^0'b\ _ f 

Hence, there must be a -i; G ker>C with f{v) 7^ as as otherwise W C ker>C and k > k. 

It remains to be shown that b ^ 7^(>C), and we will do this by contradiction: 
Suppose that b G TI{jC) with a p G X such that b = Cp. Note that for any given 
q; G M, it is also b = C{p + av)^ where v G ker £ such that f{v) ^ 0. Choose a such 
that a ^ {d- f{p))/f{v) and let p := p + av and S := {{w^^ - ^O^li'^i with 
•= — d). It can be checked that S is linearly independent by taking 

arbitrary {Pi}^^^ C M and demanding 

^ ^^(*) _ ^.p 



The second component yields = X]i=i A^i? which results in 

The set {ty^^^ljLj is, however, hnearly independent such that all /3j must vanish. Hence 
S is linearly independent. But S is also a subset of ker£ as 

This means that k > k, which is a contradiction. 
□ 
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